%n=11;
n = 101;

dt = 0.000045;
steps = 1000;

C = toeplitz(sparse([2 -1 zeros(1,n-3) -1]));

Uinit = zeros(n, steps+1);
Uinit(:,1) = (1:n)/n;

%Forward Euler

UFE = Uinit;
M = speye(n) - dt*n^2*C;

for t=2:steps+1
	UFE(:,t) = M*UFE(:,t-1);
end

%RK2

URK2 = Uinit;

for t=2:steps+1
	k = URK2(:,t-1) - dt*n*n*C*URK2(:,t-1);
	URK2(:,t) = URK2(:,t-1) - dt*n*n*C*(URK2(:,t-1) + k)/2;
end

%RK4

URK4 = Uinit;

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

%Trapezoidal

UTrap = Uinit;
M1 = eye(n) - dt*n^2*C/2;
M2 = inv(eye(n) + dt*n^2*C/2);

for t=2:steps+1
	UTrap(:,t) = M2*(M1*UTrap(:,t-1));
end

%Adams-Bashforth

UAB = Uinit;
M1 = eye(n) - 3*dt*n^2*C/2;
M2 = dt*n^2*C/2;

UAB(:,2) = (eye(n)-dt*n^2*C)*UAB(:,1);
for t=3:steps+1
	UAB(:,t) = M1*UAB(:,t-1) + M2*UAB(:,t-2);
end