function [ut]=derv(neqn,t,u)
%
% Function derv computes the derivative vector 
% for the nonlinear PDE
%
% Problem parameters
  [alpha,k,sigma,L,ui,ua]=par1;
  [ a, e, rhos,cps,ls,cs]=par2;
%
% Spatial grid
  dx=1.0/(neqn-2);
  dxs=dx*dx;
%
% Insulation
  for i=1:neqn-2
    if(i==1)
      u0=u(2)+2.0*dx*L*(sigma/k)*(a*ua^4-e*u(1)^4);
      ut(1)=(u(2)-2.0*u(1)+u0)/dxs;
    else
      ut(i)=(u(i+1)-2.0*u(i)+u(i-1))/dxs;
    end
  end
%
% Steel
  ut(neqn)=(k*L/alpha)*(1.0/cs)*(u(neqn-2)-u(neqn-1))/dx;
  ut(neqn-1)=ut(neqn);
%
% End of derv
