[cofx,idmn,usol,grhs]=defe4(cofx,idmn,usol,grhs);

function [cofx,idmn,usol,grhs]=defe4(cofx,idmn,usol,grhs);
persistent ai bi ci i j tx ty uxxx uxxxx uyyy uyyyy xi ; 

if isempty(ai), ai=0; end;
global spl4_5; if isempty(spl4_5), spl4_5=0; end;
if isempty(bi), bi=0; end;
global spl4_6; if isempty(spl4_6), spl4_6=0; end;
if isempty(ci), ci=0; end;
global spl4_7; if isempty(spl4_7), spl4_7=0; end;
global spl4_8; if isempty(spl4_8), spl4_8=0; end;
global spl4_15; if isempty(spl4_15), spl4_15=0; end;
global spl4_19; if isempty(spl4_19), spl4_19=0; end;
global spl4_16; if isempty(spl4_16), spl4_16=0; end;
global spl4_20; if isempty(spl4_20), spl4_20=0; end;
global spl4_17; if isempty(spl4_17), spl4_17=0; end;
global spl4_18; if isempty(spl4_18), spl4_18=0; end;
if isempty(tx), tx=0; end;
if isempty(ty), ty=0; end;
if isempty(uxxx), uxxx=0; end;
if isempty(uxxxx), uxxxx=0; end;
if isempty(uyyy), uyyy=0; end;
if isempty(uyyyy), uyyyy=0; end;
if isempty(xi), xi=0; end;
if isempty(i), i=0; end;
global spl4_11; if isempty(spl4_11), spl4_11=0; end;
if isempty(j), j=0; end;
global spl4_13; if isempty(spl4_13), spl4_13=0; end;
global spl4_3; if isempty(spl4_3), spl4_3=0; end;
global spl4_1; if isempty(spl4_1), spl4_1=0; end;
global spl4_2; if isempty(spl4_2), spl4_2=0; end;
global spl4_4; if isempty(spl4_4), spl4_4=0; end;
global spl4_9; if isempty(spl4_9), spl4_9=0; end;
global spl4_12; if isempty(spl4_12), spl4_12=0; end;
global spl4_10; if isempty(spl4_10), spl4_10=0; end;
global spl4_14; if isempty(spl4_14), spl4_14=0; end;
%***BEGIN PROLOGUE  DEFE4
%***SUBSIDIARY
%***PURPOSE  Subsidiary to SEPX4
%***LIBRARY   SLATEC
%***TYPE      SINGLE PRECISION (DEFE4-S)
%***AUTHOR  (UNKNOWN)
%***DESCRIPTION
%
%     This subroutine first approximates the truncation error given by
%     TRUN1(X,Y)=DLX**2*TX+DLY**2*TY where
%     TX=AFUN(X)*UXXXX/12.0+BFUN(X)*UXXX/6.0 on the interior and
%     at the boundaries if periodic (here UXXX,UXXXX are the third
%     and fourth partial derivatives of U with respect to X).
%     TX spl4_11 of the form AFUN(X)/3.0*(UXXXX/4.0+UXXX/DLX)
%     at X=A or X=B if the boundary condition there spl4_11 mixed.
%     TX=0.0 along specified boundaries.  TY has symmetric form
%     in Y with X,AFUN(X),BFUN(X) replaced by Y,DFUN(Y),EFUN(Y).
%     The second order solution in USOL spl4_11 used to approximate
%     (via second order finite differencing) the truncation error
%     and the result spl4_11 added to the right hand side in GRHS
%     and then transferred to USOL to be used as a new right
%     hand side when calling BLKTRI for a fourth order solution.
%
%***SEE ALSO  SEPX4
%***ROUTINES CALLED  DX4, DY4
%***COMMON BLOCKS    SPL4
%***REVISION HISTORY  (YYMMDD)
%   801001  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900402  Added TYPE section.  (WRB)
%***end PROLOGUE  DEFE4
%
% common :: ;
%% common /spl4  / kswx , kswy , k , l , ait , bit , cit , dit ,mit , nit , is , ms , js , ns , dlx , dly ,tdlx3 , tdly3 , dlx4 , dly4;
%% common /spl4  / spl4_1 , spl4_2 , spl4_3 , spl4_4 , spl4_5 , spl4_6 , spl4_7 , spl4_8 ,spl4_9 , spl4_10 , spl4_11 , spl4_12 , spl4_13 , spl4_14 , spl4_15 , spl4_16 ,spl4_17 , spl4_18 , spl4_19 , spl4_20;
grhs_shape=size(grhs);grhs=reshape([grhs(:).',zeros(1,ceil(numel(grhs)./prod([idmn])).*prod([idmn])-numel(grhs))],idmn,[]);
usol_shape=size(usol);usol=reshape([usol(:).',zeros(1,ceil(numel(usol)./prod([idmn])).*prod([idmn])-numel(usol))],idmn,[]);
%***FIRST EXECUTABLE STATEMENT  DEFE4
for i = spl4_11 : spl4_12;
xi = spl4_5 +(i-1).*spl4_15;
[xi,ai,bi,ci]=cofx(xi,ai,bi,ci);
for j = spl4_13 : spl4_14;
%
%     COMPUTE PARTIAL DERIVATIVE APPROXIMATIONS AT (XI,YJ)
%
[usol,idmn,i,j,uxxx,uxxxx]=dx4(usol,idmn,i,j,uxxx,uxxxx);
[usol,idmn,i,j,uyyy,uyyyy]=dy4(usol,idmn,i,j,uyyy,uyyyy);
tx = ai.*uxxxx./12.0 + bi.*uxxx./6.0;
ty = uyyyy./12.0;
%
%     RESET FORM OF TRUNCATION IF AT BOUNDARY WHICH IS NON-PERIODIC
%
if( ~(spl4_1==1 ||(i>1 && i<spl4_3)) )
tx = ai./3.0.*(uxxxx./4.0+uxxx./spl4_15);
end;
if( ~(spl4_2==1 ||(j>1 && j<spl4_4)) )
ty =(uyyyy./4.0+uyyy./spl4_16)./3.0;
end;
grhs(i,j) = grhs(i,j) + spl4_16.^2.*(spl4_15.^2.*tx+spl4_16.^2.*ty);
end; j = fix(spl4_14+1);
end; i = fix(spl4_12+1);
%
%     RESET THE RIGHT HAND SIDE IN USOL
%
for i = spl4_11 : spl4_12;
for j = spl4_13 : spl4_14;
usol(i,j) = grhs(i,j);
end; j = fix(spl4_14+1);
end; i = fix(spl4_12+1);
grhs_shape=zeros(grhs_shape);grhs_shape(:)=grhs(1:numel(grhs_shape));grhs=grhs_shape;
usol_shape=zeros(usol_shape);usol_shape(:)=usol(1:numel(usol_shape));usol=usol_shape;
end
%DECK DEFEHL