% Declare the relevant parameters
h=0.2;
w0=0.0;
tmin=0.0;
tmax=2.0;

nmax=(tmax-tmin)/h;

tlast=tmin;
wlast_euler=w0;
wlast_midpoint=w0;
wlast_rk4=w0;

t(1) = tmin;
w_euler(1) = w0;
w_midpoint(1) = w0;
w_rk4(1) = w0;
y(1) = w0;

% Beginning of Loop

for i=2:(nmax+1)

% Euler Method
    w_euler(i) = wlast_euler + h*deriv(tlast,wlast_euler);

% Midpoint Method
    a1 = wlast_midpoint + h/2*deriv(tlast,wlast_midpoint);
    w_midpoint(i) = wlast_midpoint + h*deriv(tlast + h/2,a1);

% Runge Kutta 4
    k1 = h*deriv(tlast,wlast_rk4);
    k2 = h*deriv(tlast+h/2,wlast_rk4+k1/2.0);
    k3 = h*deriv(tlast+h/2,wlast_rk4+k2/2.0);
    k4 = h*deriv(tlast+h,wlast_rk4+k3);
    w_rk4(i)=wlast_rk4+(k1 + 2*k2 + 2*k3 + k4)/6.0;

% Advance time
    t(i) = tlast + h;
    wlast_euler =w_euler(i);
    wlast_midpoint =w_midpoint(i);
    wlast_rk4 =w_rk4(i);
    tlast=t(i);

% Exact Solution
    y(i) = exact(t(i));
    d_euler(i) = abs(y(i) - w_euler(i));
    d_midpoint(i) = abs(y(i) - w_midpoint(i));
    d_rk4(i) = abs(y(i) - w_rk4(i));

end 

% Exact Solution, fine grid

  tlast=tmin;
  y_fine(1) = w0;
  h=h/10

% Beginning of Loop

  for i=2:(10*nmax+1)
    t_fine(i) = tlast + h;
    y_fine(i) = exact(t_fine(i));
    tlast=t_fine(i);
  end



% Plot Results for w, y

plot(t,w_euler,'+b');
hold on;
plot(t,w_midpoint,'og');
hold on;
plot(t,w_rk4,'xm');
hold on;
plot(t_fine,y_fine,'r');

% Label the plot 

title('Solution to dy/dt=t*exp(3*t)-2y');
xlabel('Time');
ylabel('y');
legend('Euler','Midpoint','RK4','exact')

% Save the results (for each timestep) in a file

file = zeros(size(t,2),8);

file(:,1)=t';
file(:,2)=y';
file(:,3)=w_euler';
file(:,4)=d_euler';
file(:,5)=w_midpoint';
file(:,6)=d_midpoint';
file(:,7)=w_rk4';
file(:,8)=d_rk4';

save output file -ASCII;
