% For this version, we do hybrid Runge-Kutta and Crank-Nicolson.
% Each step is much slower, but we can take much larger time steps (which 
% results in lower accuracy).

clear all;

%initial info

N = 300;

dx = 2/N;
dt = 10^-5;
x = -1:dx:1;

r1 = dt/(2*dx);
r3 = dt/(4*dx^3);
diff1 = sparse([0 r1 zeros(1,N-2) -r1]);
diff3 = sparse([0 -2*r3 r3 zeros(1,N-4) -r3 2*r3]);
d1 = toeplitz(diff1,-diff1);
d3 = 2*toeplitz(diff3,-diff3);
d3m = speye(N+1) - toeplitz(diff3,-diff3);

f1 = 200./(cosh(10.*(x+0.7)).^2);
f2 = 100./(cosh(sqrt(50)*x).^2);

U = zeros(N+1,1);
temp = U;
k1 = U;
k2 = U;
k3 = U;
k4 = U;

T = 0.015;
U = (f1+f2)';

for n=2:round(T/dt)+1
	temp = U;
	k1 = d3*(d3m\temp) + 6*temp.*(d1*temp);
	temp = U+k1/2;
	k2 = d3*(d3m\temp) + 6*temp.*(d1*temp);
	temp = U+k2/2;
	k3 = d3*(d3m\temp) + 6*temp.*(d1*temp);
	temp = U+k3;
	k4 = d3*(d3m\temp) + 6*temp.*(d1*temp);
	U = U + (k1+2*k2+2*k3+k4)/6;
end
