So, I'm writing a code to solve Hamilton's equations for a free classical particle incident on different potential barriers in 1D. I'm eventually going to be using this for some sort of semi-classical analysis, but for now have to get this working as expected. I first tested it using a Gaussian barrier, which works.. I'm now trying to incorporate a rectangular barrier, but I know I'm getting something wrong with the physics. The square barrier is just a piecewise defined as follows: U = Uo for -a<=x<=a = 0 elsewhere The Hamiltonian is as simple as it gets: H = p^2/(2m) + U So, dH/dp = p/m = dx/dt and dH/dx = dU/dx = -dp/dt These are the two equations that I propagate simultaneously using an RK4 algorithm. My problem is how to handle the square potential numerically. Since dU/dx is infinite at the edges, there's basically an infinite force acting over an infinitesimally small time dt. Any hints as to how to code this into my program? My thoughts are that it has to do with the timestep... because any arbitrarily large slope that I code in is acting over one timestep, dt, which is too large. I'm trying to think about the change in momentum analytically. I basically have the following: Δp = ∫F dt = ∫-dU/dx dt So, at the edges, dU/dx -> ∞ and dt -> 0 I know that analytically, if the incident particle does not have E > V, then it will reflect with equal and opposite momentum. So Δp = -2p = ∫-dU/dx dt I'm not quite sure how to resolve -2p in the limiting case analytically.. I think if I can figure this out, it will point me in the right direction as far as the coding goes. Any help would be greatly appreciated.