clear all;

%initial info

N = 100;
% N = 1000;
n = 2*N+1;

%Students may choose different parameters here
dt = 0.1;
% dt = 0.05;
T = 10;

u0 = [zeros(N,1);1;zeros(N,1)];
Uinit = zeros(n,T/dt+1);
Uinit(:,1) = u0;

K = toeplitz(sparse([2 -1 zeros(1,n-2)]));
I = eye(n);

%exact solution - takes a while for large N
u = expm(full(-K)*T)*u0;

%Backward Euler

UBE = Uinit;
M = I+dt*K;

for t=2:T/dt+1
	UBE(:,t) = M\UBE(:,t-1);
end

errBE = sqrt(sum((u-UBE(:,end)).^2))

%BDF2 - start with a Backward Euler step

UBDF2 = Uinit;
UBDF2(:,2) = UBE(:,2);
M = 3*I + 2*dt*K;

for t=3:T/dt+1
	UBDF2(:,t) = M\(4*UBDF2(:,t-1) - UBDF2(:,t-2));
end

errBDF2 = sqrt(sum((u-UBDF2(:,end)).^2))

% Runge-Kutta, RK2

URK2 = Uinit;

for t=2:T/dt+1
	k = URK2(:,t-1) - dt*K*URK2(:,t-1);
	URK2(:,t) = URK2(:,t-1) - dt*K*(URK2(:,t-1) + k)/2;
end

errURK2 = sqrt(sum((u-URK2(:,end)).^2))

% Runge-Kutta, RK4

URK4 = Uinit;

for t=2:T/dt+1
        k1 = -K*URK4(:,t-1);
	k2 = -K*(URK4(:,t-1) + dt*k1/2);
	k3 = -K*(URK4(:,t-1) + dt*k2/2);
	k4 = -K*(URK4(:,t-1) + dt*k3);
	URK4(:,t) = URK4(:,t-1) + dt*(k1+2*k2+2*k3+k4)/6;
end

errURK4 = sqrt(sum((u-URK4(:,end)).^2))

% ode45

[A,U45] = ode45(@Kmult,[0 T], u0);

err45 = (sqrt(sum((u-U45(end,:)').^2)))

% ode15s

[B, U15s] = ode15s(@Kmult,[0 T], u0);

err15s = (sqrt(sum((u-U15s(end,:)').^2)))