[erfresult,x]=erf(x);

function [erfresult,x]=erf(x);
erfresult=[];
persistent erf erfcs first firstCall nterf sqeps sqrtpi xbig y ; if isempty(firstCall),firstCall=1;end; 

if isempty(erfresult), erfresult=0; end;
if isempty(erfcs), erfcs=zeros(1,13); end;
if isempty(sqeps), sqeps=0; end;
if isempty(sqrtpi), sqrtpi=0; end;
if isempty(xbig), xbig=0; end;
if isempty(y), y=0; end;
if isempty(nterf), nterf=0; end;
%***BEGIN PROLOGUE  ERF
%***PURPOSE  Compute the error function.
%***LIBRARY   SLATEC (FNLIB)
%***CATEGORY  C8A, L5A1E
%***TYPE      SINGLE PRECISION (ERF-S, DERF-D)
%***KEYWORDS  ERF, ERROR FUNCTION, FNLIB, SPECIAL FUNCTIONS
%***AUTHOR  Fullerton, W., (LANL)
%***DESCRIPTION
%
% ERF(X) calculates the single precision error function for
% single precision argument X.
%
% Series for ERF        on the interval  0.          to  1.00000D+00
%                                        with weighted error   7.10E-18
%                                         log weighted error  17.15
%                               significant figures required  16.31
%                                    decimal places required  17.71
%
%***REFERENCES  (NONE)
%***ROUTINES CALLED  CSEVL, ERFC, INITS, R1MACH
%***REVISION HISTORY  (YYMMDD)
%   770401  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890531  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900727  Added EXTERNAL statement.  (WRB)
%   920618  Removed space from variable name.  (RWC, WRB)
%***end PROLOGUE  ERF
if isempty(first), first=false; end;
if firstCall,   erfcs(1)=[-.049046121234691808e0];  end;
if firstCall,   erfcs(2)=[-.14226120510371364e0];  end;
if firstCall,   erfcs(3)=[.010035582187599796e0];  end;
if firstCall,   erfcs(4)=[-.000576876469976748e0];  end;
if firstCall,   erfcs(5)=[.000027419931252196e0];  end;
if firstCall,   erfcs(6)=[-.000001104317550734e0];  end;
if firstCall,   erfcs(7)=[.000000038488755420e0];  end;
if firstCall,   erfcs(8)=[-.000000001180858253e0];  end;
if firstCall,   erfcs(9)=[.000000000032334215e0];  end;
if firstCall,   erfcs(10)=[-.000000000000799101e0];  end;
if firstCall,   erfcs(11)=[.000000000000017990e0];  end;
if firstCall,   erfcs(12)=[-.000000000000000371e0];  end;
if firstCall,   erfcs(13)=[.000000000000000007e0];  end;
if firstCall,   sqrtpi=[1.7724538509055160e0];  end;
if firstCall,   first=[true];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  ERF
if( first )
[nterf ,erfcs]=inits(erfcs,13,0.1.*r1mach(3));
xbig = sqrt(-log(sqrtpi.*r1mach(3)));
sqeps = sqrt(2.0.*r1mach(3));
end;
first = false;
%
y = abs(x);
if( y>1. )
%
% ERF(X) = 1. - ERFC(X) FOR  ABS(X) .GT. 1.
%
if( y<=xbig )
erfresult = (abs(1.0-erfc(y)).*sign(x));
end;
if( y>xbig )
erfresult = (abs(1.0).*sign(x));
end;
else;
%
% ERF(X) = 1. - ERFC(X) FOR -1. .LE. X .LE. 1.
%
if( y<=sqeps )
erfresult = 2.0.*x./sqrtpi;
end;
if( y>sqeps )
erfresult = x.*(1.0+csevl(2..*x.^2-1.,erfcs,nterf));
end;
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',x); evalin('caller',[inputname(1),'=FUntemp;']); end
return;
end;
%
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',x); evalin('caller',[inputname(1),'=FUntemp;']); end
end
%DECK EXBVP