implicit.m

clf
hold on
disp("Note: You must define the CX output as the general solution here after running one time for this code to work")
syms y(x) x d C2 C4

xrange = -5:0.1:5;
yrange = -10:0.1:10;

% Quiver
%ode = @(y, x) (x+2*y^2)/(2*y*x);
ode = @(y, x) (y+x)/(y-x);
odev = vectorize(ode);
odev = eval(odev);
[X,Y] = meshgrid(xrange, yrange);
S = feval(odev, Y, X); 
L = sqrt(1+S.^2);
quiver(X, Y, 1./L, S./L, 0.5);

% General Solution
%solutions = dsolve(diff(y)==(x+2*y^2)/(2*y*x));
solutions = dsolve(diff(y)==(y+x)/(y-x));
% After dsolve, we want y to be not be y(x)
syms y
disp(solutions);

% IVPs
ivp_problem = true;
ivp_y = [-2];
ivp_x = [-2];

% If not normal IVP use values of C
Cs = [0 1 2 3];

if ivp_problem
    sol_count = length(ivp_y);
else
    sol_count = length(Cs);
end

leg = {'Direction'};
plot_idx = 2;

last_sol = false;

for idx = 1:length(solutions)
   eq = y == solutions(idx);
   sol = solve(eq, C4);
   
   if sol == last_sol
       continue;
   end
   
   last_sol = sol;

   disp(sol)
       
   for ivp = 1:sol_count
             
       if ivp_problem
           C = eval(subs(sol, {y x}, {ivp_y(ivp), ivp_x(ivp)}));
       else
           C = Cs(ivp);
       end
       
       color = rand(1,3);
       
       disp(sol);
              
       fcontour(sol, [xrange(1) xrange(length(xrange)) yrange(1) yrange(length(yrange))],'LineWidth',2,'LineColor',color,'LevelList', [C]);
       
       if ivp_problem
           leg{plot_idx} = sprintf("%s: y(%.2f) == %.2f", char(sol), ivp_y(ivp), ivp_x(ivp));
       else
           leg{plot_idx} = sprintf("%s: C == %.2f",char(sol), C);
       end
       plot_idx = plot_idx + 1;
   end
end

legend(leg, 'FontSize', 10)
title(sprintf("ODE: %s", char(ode)));
axis([xrange(1) xrange(length(xrange)) yrange(1) yrange(length(yrange))])