function price = HopPut(S0,K,r,T,sigma,Smax,dS,dt)

M=round(Smax/dS);
dS=Smax/M;
N=round(T/dt);
dt=T/N;
matval=zeros(M+1,N+1);    
vetS=linspace(0,Smax,M+1)';
veti=0:M;
vetj=0:N;
%set up boundary conditions
matval(:,N+1) = max(K-vetS,0);
matval(1,:) = K*exp(-r*dt*(N-vetj));
matval(M+1,:) = 0;

for j=N:-1:1
    for i=2:M
        if mod(j+i,2)==1
            %Use E 
            a=0.5*dt*(sigma^2*veti-r).*veti;
            b=1-dt*(sigma^2*veti.^2+r);
            c=0.5*dt*(sigma^2*veti+r).*veti;
            matval(i,j)=a(i)*matval(i-1,j+1)+b(i)*matval(i,j+1)+c(i)*matval(i+1,j+1);
        end
    end
    for i=2:M
        if mod(j+i,2)==0
            %Use I

            x=0.5*(r*dt*veti-sigma^2*dt*(veti.^2));
                y=1+sigma^2*dt*(veti.^2)+r*dt;
                z=-0.5*(r*dt*veti+sigma^2*dt*(veti.^2));
                matval(i,j)=(1/y(i))*matval(i,j+1)-(z(i)/y(i))*matval(i+1,j)-(x(i)/y(i))*matval(i-1,j);
        end
    end
end
price=interp1(vetS,matval(:,1),S0)
