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

@author: sam
"""

import numpy as np
import matplotlib.pyplot as plt



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

def rk2solve(f,x0,y0,xmax,h): # same thing i wrote Q2's code as a function already
    x = np.array([x0])
    y = np.array([y0])
    while x0 < xmax:
        k1 = f(x0,y0)
        k2 = f(x0+(h/2),y0+(k1*h)/2)
        
        x = np.vstack([x,x0])
        y = np.vstack([y,y0])
        x0 += h
        y0 += h*k2
    return x,y



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



plt.plot(x,y)