%%%this function uses Euler's method to solve a problem posed in class
%%%on September 12. Namely solving the initial value problem
%%% x' = -y, y'=x with initial conditions
%%% x(2)=4, y(2)=0, to then figure out x(8) and y(8).
%%%
%%%The following actually allows you to pick any initial time t0 and
%%%final time tmax, any intial values x0 and y0, and any number N of steps.
%%% For the specific example above, say with N=50, you would run
%%%   Circle(2,8,50,4,0)
%%%
%Try playing around with different values of N to see what happens!
%%%
function M = Circle(t0, tmax,N,x0,y0)
Deltat = (tmax-t0)/N; %this determines the size of the time step h
x = zeros(1,N+1); %define an empty vector of x's
y = zeros(1,N+1); %define an empty vector of y's
t = zeros(1,N+1); %define an empty vector of times.
xprime = zeros(1,N);
yprime=zeros(1,N);
Deltax = zeros(1,N);
Deltay = zeros(1,N);
Newx = zeros(1,N);
Newy = zeros(1,N);
x(1) = x0; %first value in the x vector is the inital cond.
y(1) = y0; %first value in the y vector is the inital cond.
t(1)=  t0; %Our starting time is t0.
%%%
% The following loop is written in painstaking detail. If we wanted, we could
% combine steps to write x(n+1) = x(n) - Deltat*y(n), etc. 
%%% 
for n = 1:N %this is the iterative loop for Euler's Method
     xprime(n)= -y(n); % This is the equation x'=-y
     yprime(n) = x(n); % And the equation  y'=x
     Deltax(n) = Deltat*xprime(n); %How much x is changing in this step
     Deltay(n) = Deltat*yprime(n); %How much y is changing in this step
     Newx(n) = x(n) + Deltax(n); %Add the change to get the new value
     x(n+1) = Newx(n); %And record it on the next row
     Newy(n) = y(n) + Deltay(n); %Same for y
     y(n+1) = Newy(n); 
t(n+1)=t(n)+Deltat; %Finally, update the time
end %end the loop
%%%
% Now that we have our numbers, we need to plot them
%%%
h1 = figure;%provides a handle for a plot containing both x and y
%the code below is strictly used for plotting
figure(h1)
plot(t,x, 'blue', t,y, 'red')
hold on
title('Plot of both x and y')
ylabel('value'), xlabel('time')
legend('x','y')
     M = [t',x',y']; % package our data into an easy-to-read matrix.

