function asiancontinuous(S0,K,r,vol,T)

%Implementation of Vecer's PDE Method
%Vecer, J. (2002): "Unified Asian Pricing",  Risk, Vol. 15, No. 6, 113-116 


%%%%%%%%%%%%%%%%%%% MESH %%%%%%%%%%%%%%%%%%%%%%%

N = 200;       %number of subintervals in space 
M = 200;       %number of subintervals in time

%more points -> higher precision, but slower

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Xmesh  x
xmin = -1;
xmax = 1;
x = linspace(xmin, xmax, N+1);

%Tspan
t = linspace(0, T, M+1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
m = 0;
sol = pdepe(m, @pdef, @pdeic, @pdebc, x, t);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Output of the value of the option
X_0 = (1-exp(-r*T))*S0/r/T - exp(-r*T)*K;
x0 = X_0/S0;

uout = pdeval(m,x,sol(M+1,:),x0);
price = uout*S0;

fprintf( '\n\nPrice of Asian Option is %8.6f\n\n', price);

%---------------------------------------------------
%---Description of PDE---
function [c, f, s] = pdef(x, t, u, DuDx)
global vol r T;
c = 1;

    f = 0.5*vol^2*( (1-exp(-r*t))/(r*T) - x )^2*DuDx;
    s = vol^2*((1-exp(-r*t))/(r*T) - x)*DuDx;

%---Initial Condition---
function u0 = pdeic(x)
u0 = max(x, 0);

%---Boundary Condition---
function [pl, ql, pr, qr] = pdebc(xl, ul, xr, ur, t)
pl = ul;
ql = 0;
pr = ur - xr;
qr = 0;