clear;
%  Code for test problem from class example
%%%%for case where  
%%%%dy/dx=y+2x-x^2 (Greenberg equation 2, page 293
%%%%xlen is maximum x value
xlen=2.0;
%%%%dx is the grid spacing, set to=0.1
dx=0.1;
% x is the vector of grid points x(1)=0, x(2)=dx, x(nx)=xlen
x=[0:dx:xlen];
%%nx is the number of points in x
nx=length(x);
%%%% initial condition y(x=0)=1 in Greenberg.  
%%%%In matlab, we set x=0 to occur at point 1, so initial condition is y(1)=1
%% set initial condition for Euler
y(1)=1;
%set initial condition for 2nd order R-K
yrk2(1)=1;
%set initial condition for 4th order R-K
yrk4(1)=1;

%%%%%  This loop does the Euler, 2nd-order RK and 4th order RK 
for k=1:nx-1

%%%%Euler (EQN 6 page 295 Greenberg)
y(k+1)=y(k)+dx*(y(k)+2*x(k)-x(k)^2);

%%%%2nd order RK
%%xk1 computes k1, in EQN 23 Page 303 Greenberg
xk1=dx*(yrk2(k)+2*x(k)-x(k)^2);
%%xk2 computes k2 in EQN 23
xk2=dx*(yrk2(k)+xk1+2*x(k+1)-x(k+1)^2);
%%yrk2 computes y_n+1 in EQN 23
yrk2(k+1)=yrk2(k)+0.5*(xk1+xk2);

%%%%%4th order RK:  Page 305, EQN 25 Greenberg
xtemp=x(k)+dx/2;
%%compute k1
xk14=dx*(yrk4(k)+2*x(k)-x(k)^2);
%%compute k2
xk24=dx*(yrk4(k)+0.5*xk14+2*xtemp-xtemp^2);
%%compute k3
xk34=dx*(yrk4(k)+0.5*xk24+2*xtemp-xtemp^2);
%%compute k4
xk44=dx*(yrk4(k)+xk34+2*x(k+1)-x(k+1)^2);
%%%%now can compute y_n+1 in EQN 25
yrk4(k+1)=yrk4(k)+(xk14+2*xk24+2*xk34+xk44)/6;

end

%the point k+1 is the last point
k+1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  increase number of grid points/halve dx
%dx1 is the new grid spacing=0.05
dx1=0.05;
x1=[0:dx1:xlen];
nx1=length(x1)
y1(1)=1;
yrk21(1)=1;
yrk41(1)=1;
for k1=1:nx1-1

%Euler
y1(k1+1)=y1(k1)+dx1*(y1(k1)+2*x1(k1)-x1(k1)^2);

%2nd order RK
xk11=dx1*(yrk21(k1)+2*x1(k1)-x1(k1)^2);
xk21=dx1*(yrk21(k1)+xk11+2*x1(k1+1)-x1(k1+1)^2);
yrk21(k1+1)=yrk21(k1)+0.5*(xk11+xk21);

%4th order RK
xtemp1=x1(k1)+dx1/2;
xk141=dx1*(yrk41(k1)+2*x1(k1)-x1(k1)^2);
xk241=dx1*(yrk41(k1)+0.5*xk141+2*xtemp1-xtemp1^2);
xk341=dx1*(yrk41(k1)+0.5*xk241+2*xtemp1-xtemp1^2);
xk441=dx1*(yrk41(k1)+xk341+2*x1(k1+1)-x1(k1+1)^2);
yrk41(k1+1)=yrk41(k1)+(xk141+2*xk241+2*xk341+xk441)/6;

end
%the point k1+1 is the last point
k1+1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%Compute Exact Solution (which I solved in class--see lecture notes)
yexact=x1.^2+exp(x1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%PRINT RESULTS

%print all numbers on a line, change to format long (to show lots of significant figures: default is format short)
format long
aout=[x1(k1+1),yexact(nx1),y(nx),y1(nx1),yrk2(nx),yrk21(nx1),yrk4(nx),yrk41(nx1)]

%or here is a way to make nice formatted output, but be careful about the formatting issues
%(%6.4f is a number that is 7 figures long, with 4 numbers after the decimal point.  Note this will not 
%work if y >=100, then would have to change to %8.4f
sprintf('x=%3.1f yexact=%7.4f y=%7.4f y1=%7.4f yrk2=%7.4f yrk21=%7.4f yrk4=%7.4f yrk41=%7.4f', ...
x1(k1+1),yexact(nx1),y(nx),y1(nx1),yrk2(nx),yrk21(nx1),yrk4(nx),yrk41(nx1))
clf
plot(x1,yexact,'k');
hold on

plot(x,y,'*');