clear all
close all

syms x e

uexac  = e^2/2*(exp(-1/e^2) - exp(-x.^2/e^2)) + sqrt(pi)*e/2*(erf(1/e) - x.*erf(x/e));

du = diff(uexac, x)
ee = 10;
xexac = [-1:0.001:1];

uexac2  = ee^2/2*(exp(-1/ee^2) - exp(-xexac.^2/ee^2)) + sqrt(pi)*ee/2*(erf(1/ee) - xexac.*erf(xexac/ee));

duF = matlabFunction(du);

duexac = duF(xexac,ee);


figure(1)
plot(xexac, duexac,'-r')
title('Derivative')

figure(2)
plot(xexac,uexac2,'-b')
title('Function')

%du = xexac./exp(xexac.^2/ee^2) - (3991211251234741*ee*(erf(xexac./ee) + (2*xexac)./(pi^(1/2)*ee*exp(xexac.^2/ee^2))))/4503599627370496;