I have done this before. It was a combination of FORTRAN for the numbers and Tcl/Tk for the plotting. How did it work?
First, you need a trim routine which computes the forces in all three axes and interprets the remainder as an acceleration. This requires an aerodynamic model and an engine model which could be represented by tables or discrete equations, as you like. Now you will have values which feed an integration where the conditions are updated for small time steps. Rinse and repeat with a time step of 0.5 to 2 seconds.
At some point the ground speed has reached v$_{rot}$ and your numeric model has to add negative elevator deflection to lift the nose. I used a standard pitch rate of 5° per second, but maybe you want to use less. Here the time step should be reduced for improved accuracy. If lift exceeds weight, the aircraft lifts off the runway and you have the distance and time for the ground run. Now the aircraft flies and climb rate needs to be chosen such that the aircraft arrives at the obstacle height at both 1.3 times v$_{Stall}$ and the obstacle height (might be 35 or 50 ft, depending on takeoff rules). Here you get the values for the full takeoff-length and time.
How would you get the climb rate right? I simply compared the projected time until obstacle height is reached and until 1.3 times v$_{Stall}$ is reached and adjusted the flight path angle such that both became equal.
Now you need to vary all the parameters of interest and repeat the calculation for variations in:
- ambient temperature
- take-off mass
- flap setting
- runway inclination
- wind
and whatever else comes to mind. For a landing chart you reverse the sequence but the approach is pretty much the same.