function [ sol, tresid, resid, eback ] = OPExample( n )
% OPExample --- solve Ey' = -Fy + b, dimension n, and
%               return the structured backward error with 
%               the solution
%   
    D = sparse(1:n,1:n,2*ones(1,n),n,n);
    S = sparse(2:n,1:n-1,ones(1,n-1),n,n);
    E = S+D+S';
    opts = odeset( 'Mass', E, 'MStateDependence', 'none', 'MassSingular', 'no', 'RelTol', 1.0e-7, 'AbsTol', 1.0e-7 );
    function [DYP] = defun(DT,YD)
        DYP = -2*YD + [YD(2:end);0] + [0; YD(1:end-1)] + [sin(DT); zeros(n-1,1)] ;
    end

    sol = ode15s(@defun, [0,5], zeros(n,1), opts );
    [tresid, resid, yhat, yphat ] = ODEResidualMass( sol, E, @defun, 5 ); % 5 pts per step
    nump = length(tresid);
    eback = zeros( size(tresid) );
    lambda = zeros(n,1);
    for j=1:nump,
        if resid(j)~=0,
          lambda(1) = abs(resid(j))/(2*abs(yphat(1,j)) + abs(yphat(2,j)) ...
                          + 2*abs(yhat(1,j)) + abs(yhat(2,j)) + 1);
          for i=2:n-1,
              lambda(i) = abs(resid(j))/(abs(yphat(i-1,j))+2*abs(yphat(i,j)) + ...
                              abs(yphat(i+1,j)) + abs(yhat(i-1,j)) + 2*abs(yhat(i,j)) ...
                              + abs(yhat(i+1,j)) );
          end
          lambda(n) = abs(resid(j))/(2*abs(yphat(n,j)) + abs(yphat(n-1,j)) ...
                          + 2*abs(yhat(n,j)) + abs(yhat(n-1,j)) );
        end;
        eback(j) = max( lambda ) ;
    end
        
end

