%
%  math495/stat490 -- 26 nov 03 -- djm
%
%  w13diff.m:  

clear;  close all

%  walk parameters
T = 2*pi;  sig = 1.0;
Nst = 20;  dt = T/Nst;  dx = sig*dt;  t = (0:Nst)*dt;
%Nst = 500;  dt = T/Nst;  dx = sig*dt;  t = (0:Nst)*dt;

%  two cases
choice = 2;
switch choice
    case {1}
        %  initial condition
        Y0 = 0;        
        g    = cos(t);  
        mu = 0.0;  %  no decay
    case {2}
        Y0 = 0;        
        g    = ones(size(t));  
        mu = 1.0;
end

%  initialize random number generator
phone = 2914814;
rand('state',phone);  %  initialize with my office phone number

%  iterative difference equation
Ysol = [Y0 ; zeros(Nst,1)];

for k=1:Nst
    Ysol(k+1) = (1 - mu*dx)*Ysol(k) + dx*g(k+1);
end

%  plot solution
figure(1);  clf;  hold on
plot(t,Ysol,'r')
plot(t,Ysol,'k.')
axis([0 T 1.1*min(Ysol) 1.1*max(Ysol)])
ylabel('\bf Y^n in black & continuum limit Y(T) in blue')
xlabel('\bf T-axis')

if(Nst==20)
    disp('waiting ...')
    pause
end

%  plot continuum limit
switch choice
    case {1}
        plot(t,sig*sin(t),'b')
    case {2}
        plot(t,sig*(1-exp(-mu*t)),'b')
end