Finite Difference BS Call PDE

#include <iostream>
#include <fstream>
#include <cmath>
#include <vector>
using namespace std;

double explicitCallOption(double S0, double X, double T, double r, double sigma, int iMax, int jMax) {
  // declare and initialise local variables (ds,dt)
  double S_max=2*X;
  double dS=S_max/jMax;
  double dt=T/iMax;

  // create storage for the stock price and option price (old and new)
  vector<double> S(jMax+1),vOld(jMax+1),vNew(jMax+1);

  // setup and initialise the stock price
  for(int j=0; j<=jMax; j++){
    S[j] = j*dS;
  }

  // setup and initialise the final conditions on the option price
  for(int j=0; j<=jMax; j++){
    vOld[j] = max(S[j]-X,0.);
    vNew[j] = max(S[j]-X,0.);
  }

  // loop through time levels, setting the option price at each grid point, and also on the boundaries
  for(int i=iMax-1; i>=0; i--){
    // apply boundary condition at S=0
    vNew[0] = 0.;

    for(int j=1; j<=jMax-1; j++){
      double A,B,C;
      A=0.5*sigma*sigma*j*j*dt+0.5*r*j*dt;
      B=1.-sigma*sigma*j*j*dt;
      C=0.5*sigma*sigma*j*j*dt-0.5*r*j*dt;
      vNew[j] = 1./(1.+r*dt)*(A*vOld[j+1]+B*vOld[j]+C*vOld[j-1]);
    }

    // apply boundary condition at S=S_max
    vNew[jMax] = S[jMax] - X*exp(-r*(T-i*dt));

    // set old values to new
    vOld=vNew;
  }

  // get j* such that S_0 \in [ j*dS , (j*+1)dS ]
  int jstar;
  jstar = S0/dS;
  double sum=0.;
  // run 2 point Lagrange polynomial interpolation
  sum = sum + (S0 - S[jstar+1])/(S[jstar]-S[jstar+1])*vNew[jstar];
  sum = sum + (S0 - S[jstar])/(S[jstar+1]-S[jstar])*vNew[jstar+1];
  return sum;
}

int main(){
  // declare and initialise Black Scholes parameters
  double S0=1.639,X=2.,T=1.,r=0.05,sigma=0.4;

  // declare and initialise grid paramaters
  int iMax=4,jMax=4;

  cout << explicitCallOption(S0,X,T,r,sigma,iMax,jMax) << endl;
}