#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Oct 15 14:34:49 2024

@author: sam
"""

import numpy as np
import matplotlib.pyplot as plt



def f(x,y):
    return x*(1+4*y**2)

def rk(f,x0,y0,xmax,h): # 2nd order Runge-Kutta method
    x = np.array([x0])
    y = np.array([y0])
    while x0 < xmax:
        k1 = f(x0,y0) # estimates midpoint using Euler step
        k2 = f(x0+(h/2),y0+(k1*h)/2) # computes slope at estimated point
        
        x = np.vstack([x,x0]) # similar stuff
        y = np.vstack([y,y0])
        x0 += h
        y0 += h*k2
    return x,y



x,y = rk(f,0,0,1,0.1)



plt.plot(x,y)