% Define the boundaries and system parameters
XMin = -5;
XMax = 5;

% Has to do with viewing rectangle
YMin = -.5;
YMax = .5;
YDotMin = -.5;
YDotMax = .5;

% Duration of simulation
TMax = 4;

v = 1;

% Numerical parameters
dt = 0.025;
dx = 0.01;

% Note the relationship between x and Matlab array index:
% x = XMin + (index-1)*dx
N = floor( (XMax - XMin) / dx );
x = linspace( XMin, XMax, N );
t = [0:dt:TMax];

% Initialize dynamic variables;
y = zeros(1,N);

dydx = zeros(1,N);
dydx2 = zeros(1,N);

dydt = zeros(1,N);
dydt2 = zeros(1,N);

% Initial conditions;
x_0 = 1;
a = 0.8;

for i = 1:N
    if (x(i) > 0) & (x(i) < x_0)
        dydt(i) = a*[(x(i)/x_0) - (x(i)/x_0)^5];
    end
end

% Dynamics begins here;
for j = 1:length(t)
    % Plot the current state
    subplot(211);
    plot(x,dydt);
    title('Transverse velocity: dy/dt');
    xlabel('x');
    ylabel('dy/dt');
    axis([XMin XMax YDotMin YDotMax]);
    subplot(212);
    plot(x,y);
    title('Displacement of string: y(x)');
    xlabel('x');
    ylabel('y');
    axis([XMin XMax YMin YMax]);
    
    hold on;
    
    % Update for next frame;
    % The string
    % smearing parameter s
    s = 1;
    for i = (1+s):N-s
       y(i) = y(i) + 1/2*(dydt(i-s)+dydt(i+s))*dt;
    end
    
    % Compute dydx
    % smearing parameter s
    s = 1;
    for i = (1+s):N-s
      dydx(i) = [y(i+s) - y(i-s)]/((2*s+1)*dx);
    end
    
    %plot(x,dydx,'r');
    
    % Compute dydx2
    s = 2;
    for i = (1+s):N-s
      dydx2(i) = [dydx(i+s) - dydx(i-s)]/((2*s+1)*dx);
    end
    
    %plot(x,dydx2,'g');
    axis([XMin XMax YMin YMax]);
    hold off;
    drawnow;
    
    % Relate dydx2 to dydt2 via wave equation
    dydt2 = v^2 * dydx2;
    
    % Update dydt
    for i = 1:N-1
       dydt(i) = dydt(i) + dydt2(i)*dt;
    end
   
end
