delegate double RealFunction(double arg); 
   
  public double h = 10e-6;  
   
  static double Derivative(RealFunction f, double arg) 
  { 
      double h2 = h*2; 
      return (f(x-h2) - 8*f(x-h) + 8*f(x+h) - f(x+h2)) / (h2*6); 
  } 
   
   
  

 
   
   
  interface IFunction 
  { 
      // Since operator () can't be overloaded, we'll use this trick. 
      double this[double arg] { get; } 
  } 
   
  class Function : IFunction 
  { 
      RealFunction func; 
   
      public Function(RealFunction func) 
      { this.func = func; } 
   
      public double this[double arg] 
      { get { return func(arg); } } 
  } 
   
  class Derivative : IFunction 
  { 
      IFunction func; 
      public static double h = 10e-6; 
   
      public Derivative(IFunction func) 
      { this.func = func; } 
   
      public double this[double arg] 
      { 
          get 
          { 
              double h2 = h*2; 
              return ( 
                  func[arg - h2] - func[arg + h2] + 
                  ( func[arg + h]  - func[arg - h] ) * 8 
                  ) / (h2 * 6); 
          } 
      } 
  } 
  

 
  using System; 
  using System.Collections.Generic; 
  using System.Text; 
   
  namespace IntegralCalculus 
  { 
      class Integration 
      { 
   
          // What are BMR, MAP, k1? Store their values here so they can be accessed in all your functions 
          double k1 = 0; 
          string BMR = "foo"; 
          // etc. 
   
          public double Function(double t) 
          { 
              double result = /* BMR + (MAP –BMR)(1-e-t/k1) */; 
              return result; 
          } 
   
          public double Integrate(double start, double end, long intervals) 
          { 
              // Loop through values of t to call our function 
              double value = 0; 
              // Note: it would be more optimal to store the *increment* value here, and reuse it at every step of the loop 
              double increment = (end - start) / (double)intervals; 
              for (long n = 0; n < intervals; n++) 
              { 
                  // Calculate t, interpolating linearly between start and end values 
                  // Note: it would be more computationally optimal to reuse the increment value we have already stored, 
                  // however this could be less *accurate* mathematically than recalculating it every time, because on very 
                  // large values of intervals this method could improve precision ... at least, I think! 
                  // Also notice the use of (casting) needed to allow us to combine double precision values with integer ones 
                  double t = start + ((end - start) * (double)n) / (double)intervals; 
                  // Now we call the function and fetch the result 
                  double result = Function(); 
                  // Multiply the result by the increment and add it to our value 
                  value += result * increment; 
              } 
              // Return the value 
          } 
   
          public double GetPT(double tStart, double tEnd, long intervals) { 
              double finalResult = /* (A/T)(1 – e-T/k2) + Integrate(tStart,tEnd,intervals) */; 
              return finalResult; 
          } 
   
          public void Run() 
          { 
              // Run the entire function for t[0,1] at 1000000000 precision and output to console 
              Console.WriteLine(0,1,GetPT(1000000000).ToString()); 
          } 
   
      } 
  }