R code:

#The characteristic function of Black-Scholes model

cf=function(S_0,mu,sigma,T,u){

im=complex(real=0,imaginary=1)

x=exp(im*u*(log(S_0)+(mu-((sigma**2)/2))*T)-(((sigma**2)*T*(u**2))/2))

x

}

#The psi function (eq(4) in link)

psim=function(r,S_0,mu,sigma,T,alpha,v){

im=complex(real=0,imaginary=1)

u=v-(alpha+1)*im

x=(exp(-r*T)*cf(S_0,mu,sigma,T,u))/(alpha**2+alpha-v**2+im*(2*alpha+1)*v)

x

}

#function used in the simpson weights

kdf=function(x){

if(x==0){

result=1

}

else{

result=0

}

result

}

#function used in simpson weights

kdfn=function(x,N){

if(x==N){

result=3

}

else{

result=0

}

result

}

#FFT of the eq(6) in the link

callcm=function(r,S_0,mu,sigma,T,alpha,N,h){

im=complex(real=0,imaginary=1)

v=seq(0,((N-1)*h),by=h)

eta=h

b=pi/eta

lambda=(2*pi)/(N*eta)

u=seq(1,N,by=1)

k_u=-b+lambda*(u-1)

f=function(i){

sum=0

for(j in 1:N){

sum=sum+((exp((-im*lambda*eta*(j-1)*(i-1))+im*v[j]*b)*

psim(r,S_0,mu,sigma,T,alpha,v[j]))*(eta/3)*(3+(-1)**(j)-kdf(j-1)-kdfn(j,N)))

}

sum

}

result=sapply(u,f)

result=((exp(-alpha*k_u))/(pi))*Re(result)

result

}

I've been using parameter values:

N=2^11

h=0.01

S_0=210.59

T=4/365

r=0.002175

alpha=1

mu=r

sigma=0.1404

But have at best been getting results that are correct to 2 dp, whereas I am led to believe I should be getting results as accurate as say 10^-14. So any help on where you can see errors in my formulation or the code would be much appreciated thank you!