% script: example3.m
% use mu(t)=  t^2 for this example
% also use F(x,t)=0, alpha(t)=0 (the rhs of right bdy cnd),
% phi(x)=2-x/2+x^2*sin(3*x)/3 and c(x)=1+cos(3*x)/2.
% Check files alpha.m, F.m and mu.m.  Other functions
% are hardwired into this example.
% Choices of example 3 are malpha=0.2.

global mindex;
global malpha;
% files work for Octave or Matlab

el = 1 % right hand spatial limit
T = 3 % time limit
gmma = 1.0 % gamma parameter as above
bta = 1.0 % beta boundary parameter as above
M = 100 % space node indices from 0..M
N = 300 % time node indices from 0..N
dx = el/M;
xnodes = linspace(dx,el,M)';
pnodes = 1+cos(3*xnodes)/2;

% Here is an operator experiment:
mindex = 0;
dx = el/M;
x = linspace(dx,el,M)';
% this is our "measurement"
gnodes  = [mu(T);mu(T)*exp(-xnodes)];
% now construct a first approximation to q(x), namely h(x)
mindex = 1; % now mu calculates mu'(t)
mumax = 0;
tmpt = linspace(0,T,N);
for jj = 1:N
  mumax = max(mumax,abs(mu(tmpt(jj))));
end;
qnodes = ((1+pnodes(1:M-1)*dx/2).*gnodes(1:M-1) - ...
         2*gnodes(2:M) + (1-pnodes(1:M-1)*dx/2).*gnodes(3:M+1));
% this calculation assumes d=1, else stick it in as factor
qnodes = [qnodes;0]/(mumax/gmma*(1-exp(-gmma*T))*dx*dx);
% do the last entry by quadratic interpolation
qnodes(M) = 3*qnodes(M-1)-3*qnodes(M-2)+qnodes(M-3);
%clearplot
%closeplot
plot(xnodes,qnodes);
disp('qnodes(5):')
disp(qnodes(5))
hold on
for jj = 1:8
  qtmp = operatorA(M,N,el,T,gmma,bta,gnodes,pnodes,qnodes);
  norm(qtmp-qnodes,inf)
  qnodes = qtmp;
  if (jj < 5)
    plot(x,qnodes);
    disp('qnodes(5):')
    disp(qnodes(5))
  end;
end;
% what does this computed q(x) result in? Let's see...
mindex = 0; % mu now calculates mu(t)
newgnodes = [mu(T);forwardsolv(M,N,el,T,gmma,bta,pnodes,qnodes)];
disp('Here is the error in original g minus implied g:');
disp(norm(gnodes-newgnodes,inf));