rungekutta.r

# vanilla RK4 written in R (09Dec2011) v1.1
#	evaluates a 2nd order ODE by splitting into 2- 1st order ODE 
#	and solved by Runge-Kutta Method
# Concepts and algorithm from:
#  http://en.wikipedia.org/wiki/Runge-Kutta_methods
#  http://www.cms.livjm.ac.uk/etchells/notes/ma200/2ndodes.pdf
#
#	f1(t,y) : substituted function from dz(y')/dt = F(z,y,t), z = y'
#	
#	rk4(f, h, tinit, yinit, limit): main RK4 routine
#		f - function defining the differential equation of
#		    the form dy/dt = f(t,y)
#		h - interval step
#		t/y/zinit - boundary values, where y(tinit) = yinit
#		limit - number of iterations the algorithm is done
#		 
#	TODO: fix documentation, instructions; verify using test cases
#         and compare with plots using WolframAlpha

f2 <- function(z,t,y) {
	value <- z
	value
	}
f1 <- function(z,t,y) {
	value <- -t*y
	value
	}

rk4 <- function(f1,f2 , h, tinit, yinit, zinit, limit) {
	k <- c(0,0,0,0)
	j <- c(0,0,0,0)
	y <- yinit
	t <- tinit
	z <- zinit
	tplot <- c()
	yplot <- c()
	zplot <- c()
	while (limit > 0) {
		
		
		k[1] <- h*f2(z,t,y)
		k[2] <- h*f2(z,t + 0.5*h, y + 0.5*k[1])
		k[3] <- h*f2(z,t + 0.5*h, y + 0.5*k[2])
		k[4] <- h*f2(z,t + h, y + k[3])
		
		y = y + (1./6)*(k[1] + 2*k[2] + 2*k[3] + k[4])
		
		j[1] <- h*f1(z,t,y)
		j[2] <- h*f1(z,t + 0.5*h, y + 0.5*j[1])
		j[3] <- h*f1(z,t + 0.5*h, y + 0.5*j[2])
		j[4] <- h*f1(z,t + h, y + j[3])
		
		z = z + (1./6)*(j[1] + 2*j[2] + 2*j[3] + j[4])
		
		t = t + h
		
		limit <- limit - 1
		tplot[limit] = t
		yplot[limit] = y
		zplot[limit] = z
		}
	#plot(tplot,yplot)  ## plots y(t)  vs. t
	plot(yplot,zplot) ## plots y'(t) vs y(t)
		
	}

Hi, I tried to run the code in Rstudio, but it did not create a plot (opposed to other programmings I ran so far in parallel). Any ideas?

sorry for this basic question, I am a beginner in R....