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

@author: sam
"""

import numpy as np
import matplotlib.pyplot as plt



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

def rk4solve(f,x0,y0,xmax,h): # 4th order Runge-Kutta
    x = np.array([x0])
    y = np.array([y0])
    while x0 < xmax:
        k1 = f(x0,y0)
        k2 = f(x0+(h/2),y0+(k1*h)/2)
        k3 = f(x0+(h/2),y0+(k2*h)/2)
        k4 = f(x0+h,y0+k3)
        
        x = np.vstack([x,x0])
        y = np.vstack([y,y0])
        x0 += h
        y0 += (h/6)*(k1+(2*k2)+(2*k3)+k4) # increases each y step by a weighted average slop
    return x,y



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



plt.plot(x,y)