Apologies upfront for such a large question. Attached is a write up of the problem in case any further info is needed.
These are the equations I am trying to solve: $$\frac{d}{dx}U_n(x)=1-U_n(x)\frac{(r-D+a)}{a}-U^2_n(x)\frac{1}{a\Delta\tau}, \>\>\>\>\> U_n(0)=0,$$ $$\frac{d}{dx}u_n(x)=-U_n(x)\frac{1}{a\Delta\tau}u_n(x)+U_n(x)\frac{v_{n-1}(x)}{a\Delta\tau}, \>\>\>\>\> u_n(0)=1,$$ $$\frac{d}{dx}z_n(x)=\Bigg(\frac{U_n(x)}{a\Delta\tau}+\frac{(r-D+a)}{a}\Bigg)z_n(x)+\frac{(u_n(x)-v_{n-1}(x))}{a\Delta\tau}, \\ z_n(L+M)=-\frac{u_n(L+M)}{U_n(L+M)}.$$ With the end goal of finding $v_n(x)=U_n(x)z_n(x)+u_n(x)$. Initial condition for $v_n(x)$ is $v_0(x)=\text{max}\{1-e^{x-L},0\}$, $1<<L\in\mathbb{R}$. (Note $L,M$ are constants where we assume $v_n(0)=1$ and $v_n(L+M)=0$). The subscript $n$ refers to the $n$ time steps we choose for the interval $(0,t_*)$, each step has size $\Delta\tau$. Essentially we use Euler's method to discretise $t$ turning our original PDE into this system of ODEs and then solve this system at each time step, giving us an approximation of $v(x,t)$. I've heard this being reffered to as using a grid or lattice to solve a PDE.
(This system is actually taken from the paper written by Jiang, Chen, Wang and Zhang: A new well-posed algorithm to recover implied local volatility. Appendix A.1. though I've simplified them to make $a(y)=a\in\mathbb{R}$.)
I have written some code for how I would go about solving them but the results are quite wild which makes me think I've made a mistake. I think I am right in how I have solved the first two equations, but the final one needs to be integrated backwards and this got me stuck and I think it could be what's throwing me off. This is the reverse equation I have used, is it correct? Let $1<<L,M\in\mathbb{R}$ and $w_n(s)=z_n(L+M-s)$ then $$\frac{d}{dx}w_n(s)=-\Bigg(\frac{U_n(L+M-s)}{a\Delta\tau}+\frac{(r-D+a)}{a}\Bigg)w_n(s)+\frac{(u_n(L+M-s)-v_{n-1}(L+M-s))}{a\Delta\tau}, \\ w_n(0)=-\frac{u_n(L+M)}{U_n(L+M)}.$$
Finally, the code I've written in R.
tstar=1
n<-100 #determines time intervals
dt<-tstar/n
t=seq(from=0,to=tstar,by=dt)
L=100
M=100
#L,M,n and truncno (below) can all be varied to see how the plot below changes.
x=seq(from=1,to=L+M,by=0.1)
v0<-ifelse(1-exp(x-L)<0,0,1-exp(x-L)) #initial condition
#chosen constants for the problem
a=0.25
r=0.12
D=0
#blank matrix to store each successive output
resultsmatrix=matrix(data=0,nrow=length(x),ncol=length(t))
#first column is initial condition
resultsmatrix[,1]=v0
#loop for solving the system of ODEs at each time step, note t[-1] is used here since column one of resultsmatrix is already filled with the initial condition.
for(val in (1:length(t[-1]))){
p<-c(r,D)
riccati<-function(times,y,parms){
dU.dx<-1-y[1]*((p[1]-p[2]+a)/a)-(y[1]^2)/(a*dt)
return(list(dU.dx))
}
y0<-0
sol<-ode(y=y0,times=x,func=riccati,parms=p)
linear<-function(times,y,parms){
du.dx<-(-sol[[times,2]]*y[1]/(a*dt))+sol[[times,2]]*resultsmatrix[times,val]/(a*dt)
return(list(du.dx))
}
y1<-1
sol2<-ode(y=y1,times=x,func=linear,parms=NULL)
linearzreverse<-function(times,y,parms){
tt=tail(x,n=1)+1
dw.dx<-(-(sol[[tt-times,2]]/(a*dt)+(p[1]-p[2]+a)/a))*y[1]+(sol2[[tt-times,2]]-resultsmatrix[tt-times,val])/(a*dt)
return(list(dw.dx))
}
y2<-(-tail(sol2[,2],n=1)/tail(sol[,2],n=1))
sol3<-ode(y=y2,times=x,func=linearzreverse,parms=p)
resultsmatrix[,val+1]=sol[,2]*sol3[,2]+sol2[,2]
}
#Where you see tt-times, this is my attempt at dealing with the reversal of the equation.
#Columns of resultsmatrix are the solutions at each time step. Doing the following problem specific scalings gives us the following code
S=50
yvals=x-L
K=S*exp(yvals)
truncno=100
persp(x=head(K,n=truncno),y=t,z=head(S*resultsmatrix,n=truncno),
xlab="K",ylab="t",zlab="V(K,t)", theta = 50, phi = 25,
ticktype="detailed",axes=T)
#This plot is what looks to wild to me.
