function [x,t,u,error] = Test_MOL_Euler(index_u,t_max,h,m)

% Try this test program with 
%    [x,t,u,error] = Test_MOL_Euler(2,5,1/128,8);

[x,t,u] = MOL_Euler(@d0,@d1,@f,@G,t_max,h,m);

% Graph numerical solution
[X,T] = meshgrid(x,t);
figure; mesh(X,T,u); shading interp
xlabel('x'); ylabel('t');
title(['Numerical solution u: index of u = ',num2str(index_u)])
disp('Press any key to continue.'); pause 

% Graph error in numerical solution
true_u = true_soln(X,T); error = true_u - u;
disp(['Maximum error = ',num2str(max(max(abs(error))))])
figure; mesh(X,T,error); shading interp
xlabel('x'); ylabel('t');
title(['Error in numerical solution u: index of u = ',num2str(index_u)])
disp('Press any key to continue.'); pause

% Produce maximum errors over x as t varies.
maxerr_in_x = max(abs(error'));
figure; plot(t,maxerr_in_x); text(1.02*t_max,0,'t')
title('Maximum error for x in [0,1], as a function of t')

function true_u = true_soln(z,s)

switch index_u
    case 1
        true_u = s.^2 + z.^4;
    case 2
        true_u = exp(-0.1*s).*sin(pi*z);
    case 3
        true_u = exp(1./(1+s)).*cos(pi*z);
    case 4
        true_u = (1-exp(-s)).*cos(pi*z);
end
end % true_u

function answer = G(z,s)
    
% This routine assumes s is a scalar, while z can be a vector.
switch index_u
    case 1
        answer = 2*s - 12*z.^2;
    case 2
        answer = (pi^2 - 0.1)*exp(-0.1*s).*sin(pi*z);
    case 3
        answer = exp(1./(1+s)).*cos(pi*z).*(pi^2 - 1./(1+s).^2);
    case 4
        answer = cos(pi*z).*(pi^2 + (1- pi^2)*exp(-s));
end
end % G

function answer = d0(s)    
z = zeros(size(s));
answer = true_soln(z,s);
end % d0
    
function answer = d1(s)
z = ones(size(s));
answer = true_soln(z,s);
end % d1

function answer = f(z)    
s = zeros(size(z));
answer = true_soln(z,s);
end % f

end % Test_MOL_Euler