% slopefield - MECH 221 script supporting Computer Lab 04

% (c) 2014-10-23, Philip D Loewen

% Setup zone
Ndots = 14; % Bigger numbers give better pictures, but take longer

% Automatic operation zone
%  A general ODE has the form  y'=F(x,y).
%  Here is our particular function F:
F = @(x,y)  (3*x.*x + 2*x) ./ (3*y.*y - 1);

% Make and plot a grid of evaluation points, size Ndots by Ndots
[X,Y] = meshgrid(linspace(-2,2,Ndots),linspace(-2,2,Ndots));
plot(X(:),Y(:),'k.');

% Fine-tune the plot. Students should look closely at this and learn:
axis equal;
hold on;
set(gca,'XLim',[-2,2],'YLim',[-2,2]);
grid on;
xlabel('x');
ylabel('y');
title('Slope field for (3y^2-1)y'' = 3x^2 + 2x');

% Calculate the horizontal separation of adjacent grid points:
[M,N] = size(X);
dx = X(1,2)-X(1,1);

% Nested loops visit each point on the grid and draw little segments:
for i=1:M
  for j=1:N
    x0 = X(i,j);          % x-coord of this grid point
    y0 = Y(i,j);          % y-coord of this grid point
    th = atan(F(x0,y0));  % theta, from tan(theta)=y'
    x  = x0 + (dx/3)*cos(th)*[-1,1];  % x-coords of segment endpts
    y  = y0 + (dx/3)*sin(th)*[-1,1];  % y-coords of segment endpts
    plot(x,y);            % draw a short segment with correct slope
  end;
end;

% Contours get generated here:
C = 1;
g = @(x,y)  x.*x + y.*y;
contour(X,Y,g(X,Y),[C,C],'r','LineWidth',2);  % One contour, g(x,y)=C
% To add 11 more contours, auto-coloured, un-comment the next line:
% contour(X,Y,g(X,Y),11);
