C++ program has stopped working- Solving ordinary differential equations

Question

I'm writing a C++ program to find solutions for first order differential equations for a college assignment. The program starts up and then once I enter the number of iterations to do I get the error message "Euler's method.exe has stopped working". This is my code:

#include <functional>
#include <vector>

using namespace std;

 
double f_r(double x, double r) {  
  return r;
  } 

double f_s(double x, double s) {
  return -x/s;  
  }
 

double eulerstep(const function<double(double,double)>& f, double xsub0, double ysub0, double h) {
  double ysub1 = ysub0+ h * f(xsub0,ysub0);
  return ysub1;
  }


double euler(const function<double(double,double)>& f, double xsub0, double ysub0, double h, int n) {
  vector<double> xsub;
  vector<double> ysub;
  xsub[0] = xsub0;
  ysub[0] = ysub0;
  n = ysub.size();
  for (int i=1; i<n; i++){
  	xsub[i+1] = xsub[i] + h;
  	ysub[i+1] = ysub[i] + h * f(xsub[i],ysub[i]);
  	cout << xsub[i] << " , " << ysub[i] << endl;
  }
  return ysub[n];
  }
  
int main() {
  int nsteps = 0;
  cout << "Number of steps?" << endl;
  cin >> nsteps;
  double h = 1.0 / nsteps;

  double r = euler(f_r,0,1,h,nsteps);
  

  double s = euler(f_s,0,1,h,nsteps);
 
  return 0;
}

I'm fairly sure the problem lies in how I've defined my vectors but I'm new to using them so can't see where I've gone wrong. I'd be extremely grateful if anyone could point out the mistakes in my method

Thanks

Answer 1

Here's an updated function. You have to set the sizes, not ask an empty vector for its size.

double euler(const function<double(double, double)>& f, double xsub0, double ysub0, double h, int n)
{
    vector<double> xsub;
    vector<double> ysub;
    xsub.resize(n+1); // so we can access [n], it must be size n+1
    ysub.resize(n+1);
    xsub[0] = xsub0;
    ysub[0] = ysub0;
    for (int i = 1; i<n; i++) {
        xsub[i + 1] = xsub[i] + h;
        ysub[i + 1] = ysub[i] + h * f(xsub[i], ysub[i]);
        cout << xsub[i] << " , " << ysub[i] << endl;
    }
    return ysub[n];
}