This script will allow you to approximate a solution to a differential equation using Euler’s method.

When inputting your dy/dx you must explicitly define every operation

Examples of inputting different dy/dx are below:

if dy/dx = 2x/y then input 2*x/y

if dy/dx = 0.5x+y then input 0.5*x+y

Watch a tutorial at : https://youtu.be/Q9BwmdPH9-M?si=O26JjlLRm4A3En69

from matplotlib.pyplot import *
from math import *

print("Euler's Method with Slope Field")
print("Enter dy/dx explicitly")
print("Examples: 2*x/y, sqrt(x)*y, sin(x) - y")

def get_float(prompt):
    while True:
        raw = input(prompt).strip()
        try:
            return float(raw)
        except:
            print("Please enter a valid number (use decimal point).")

# allowed math functions for eval
allowed_math = {
    "sqrt": sqrt, "sin": sin, "cos": cos, "tan": tan,
    "log": log, "exp": exp, "pi": pi, "e": e,
    "asin": asin, "acos": acos, "atan": atan,
    "abs": abs, "pow": pow
}

# inputs
dydx_expression = input("dy/dx = ")
x0 = get_float("Initial x = ")
y0 = get_float("Initial y = ")
h_raw = get_float("Step size h = ")
x_end = get_float("End x = ")

# handle zero step and direction
if h_raw == 0:
    print("Step size h cannot be 0. Re-run and enter a nonzero h.")
else:
    direction = 1 if x_end > x0 else -1
    h = abs(h_raw) * direction

    def dydx(_x, _y):
        env = {"x": _x, "y": _y, "__builtins__": {}}
        env.update(allowed_math)
        try:
            return eval(dydx_expression, env)
        except Exception as e:
            print("Error in dydx:", e)
            return 0

    # Euler integration
    x_vals = [x0]
    y_vals = [y0]

    while True:
        cur_x = x_vals[-1]
        if (direction == 1 and cur_x >= x_end) or (direction == -1 and cur_x <= x_end):
            break
        step = h
        if (direction == 1 and cur_x + step > x_end) or (direction == -1 and cur_x + step < x_end):
            step = x_end - cur_x
        slope = dydx(cur_x, y_vals[-1])
        new_x = cur_x + step
        new_y = y_vals[-1] + step * slope
        x_vals.append(new_x)
        y_vals.append(new_y)
        if new_x == x_end:
            break

    # print table
    print("\nTable of values:")
    print("Step\t   x\t\t   y")
    for i in range(len(x_vals)):
        print(i, "\t", "%.5f" % x_vals[i], "\t", "%.5f" % y_vals[i])

    # wait for user before plotting
    input("\nPress Enter to show the plot...")

    # Create slope field
    x_min = int(min(x_vals) - 2*abs(h))
    x_max = int(max(x_vals) + 2*abs(h))
    y_min = int(min(y_vals) - 1)
    y_max = int(max(y_vals) + 1)

    x_points = range(x_min, x_max + 1)
    y_points = range(y_min, y_max + 1)

    d = 0.25
    xs, ys, dxs, dys = [], [], [], []

    for X in x_points:
        for Y in y_points:
            m = dydx(X, Y)
            dx = d / sqrt(1 + m**2)
            dy = m * dx
            xs.append(X - dx)
            ys.append(Y - dy)
            dxs.append(2 * dx)
            dys.append(2 * dy)

    # Plot slope field
    axis((x_min, x_max, y_min, y_max))
    for i in range(len(xs)):
        arrow(xs[i], ys[i], dxs[i], dys[i], color="red", head_width=0.05, length_includes_head=True)

    # Plot Euler curve
    plot(x_vals, y_vals)
    
    show()