This is done using Student:-NumericalAnalysis:-InitialValueProblem as follows.
Note this is called heun method in Maple.
ode:=diff(y(x),x)=2*x*y(x)+x; IC:=y(0)=1; ending_x:=2; Student:-NumericalAnalysis:-InitialValueProblem(ode,IC, x=ending_x,'method'='rungekutta','submethod'='heun', 'numsteps'=41,'output'='plot','plotoptions'=['gridlines']);
To obtain table of values instead do this
ode:=diff(y(x),x)=2*x*y(x)+x; IC:=y(0)=1; ending_x:=2; Student:-NumericalAnalysis:-InitialValueProblem(ode,IC, x=ending_x,'method'='rungekutta','submethod'='heun', 'numsteps'=41,'output'='information');