%------------------------------------------------------------
% Script name: e_integ
% Date:        3/27/2007
% Description: Attempt to integrate d2x/dt2 = k^2*x
%------------------------------------------------------------
function [x0 d_x0 error] = e_integ

N   = 20;
lim = 2;
[x0 d_x0] = meshgrid([0:lim/(N-1):lim],[-lim:lim/(N-1):0]);
error     = zeros(N,N);

for k1 = 1:N       % Fix initial derivative
    for k2 = 1:N   % The initial position
       error(k1,k2) = solve(x0(1,k2), d_x0(k1,1));
    end
end
surf(x0,d_x0,error);

% RK4 integration of d2x/dt2 = k^2*x
%------------------------------------------------------------
function error = solve(x0, d_x0)

N = 500;              % Number of steps
t = linspace(0,1,N);  % The region of integration
x = zeros(N,2);       % Solution vector

x(1,1) = x0;
x(1,2) = d_x0;

k  = 2;
dt = t(2)-t(1);

for j = 2:N
    x0 = x(j-1,:);
    
    k1 = dt * [x0(2) k^2*x0(1)^2];
    x1 = x0 + 1/2*k1;
    
    k2 = dt * [x1(2) k^2*x1(1)^2];
    x2 = x0 + 1/2*k2;
    
    k3 = dt * [x2(2) k^2*x2(1)^2];
    x3 = x0 + k3;
    
    k4 = dt * [x3(2) k^2*x3(1)^2];
    
    x(j,:) = x0 + 1/6*(k1 + 2*k2 + 2*k3 + k4);
end

% Results
%------------------------------------------------------------
% plot(t,x(:,1));
% hold on;
error = x(length(x(:,1)),1);
