I want sure rasp can generate a lot of 40khz square wave at the same time accurately
It's not clear what you mean by "a lot" here: As in, the number of consecutive waves, or the pins that can generate waves simultaneously?
If the latter, basic PWM generation on the Pi only has two independent channels, which can be accessed via 4 pins; GPIOs 12 and 18 for channel 0, 13 and 19 for channel 1.
I believe pigpio has a number of other ways of generating pulses by making use of some of the other clocks on the Pi and joan's links should lead you to examples of that. With regard to the basic hardware mentioned in the last paragraph, here's a simple C program using libbcm2835 to generate a 40 khz pulse. This should work on all models.1 Part of my point is whatever implementation simulink uses, it probably only makes use of these two basic channels.
Note most of the code is actually OS signal handling so you can start it, let it run in the foreground, then stop it with ctrl-c and have the pin reset as an input. Actually turning on the wave in hardware only take a few lines and the process then just sleeps indefinitely as it does not have to do anything except wait for SIGTERM.
#include <bcm2835.h>
#include <stdio.h>
#include <signal.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>
// Compile --std=gnu99 for sigaction.
// 19200000 / 10 / 48 = 40000 Hz
#define CLOCK_DIVISOR 10
#define PWM_RANGE 48
#define PWM_PIN 18
// For pins 13 or 19, use channel 1.
#define PWM_CHANNEL 0
// For pins 12 or 13, use ALT0.
#define PWM_ALT BCM2835_GPIO_FSEL_ALT5
void sigHandler (int sig)
{
fprintf(stderr, "Stopped by signal %d.\n", sig);
bcm2835_pwm_set_mode(PWM_CHANNEL, 0, 0);
bcm2835_gpio_clr(PWM_PIN);
bcm2835_gpio_fsel(PWM_PIN, BCM2835_GPIO_FSEL_INPT);
bcm2835_close();
exit(0);
}
void setSigHandler ()
{
int sigs[] = { SIGINT, SIGTERM, SIGQUIT, SIGHUP, SIGABRT },
num_sigs = sizeof(sigs) / sizeof(int);
struct sigaction sigact;
memset(&sigact, 0, sizeof(sigact));
sigact.sa_handler = sigHandler;
for (int i = 0; i < num_sigs; i++) {
if (sigaction(sigs[i], &sigact, NULL) == -1) {
fprintf(stderr, "sigaction(%d) fail:\n", sigs[i]);
perror(NULL);
}
}
}
int main (void)
{
if (!bcm2835_init()) {
fprintf(stderr, "Failed to initialize.\n");
return 1;
}
setSigHandler();
bcm2835_pwm_set_clock(CLOCK_DIVISOR);
bcm2835_pwm_set_mode(PWM_CHANNEL, 1, 1);
bcm2835_pwm_set_range(PWM_CHANNEL, PWM_RANGE);
bcm2835_gpio_fsel(PWM_PIN, PWM_ALT);
bcm2835_pwm_set_data(PWM_CHANNEL, PWM_RANGE / 2);
while (1) sleep (60);
return 0;
}
So to compile that:
gcc --std=gnu99 square_wave.c -lbcm2835
Here's 1 ms of output sampled with a logic analyzer working at 12 Mhz:

Looks like a perfect 40 kHz square wave. However, again, you can only generate 2 of these independently using the basic hardware, which you can access via 4 pins on the 40-pin models (there's only 1 explicit PWM pin on the older 26-pin models).
1. When I first posted this I said it did not work on the Pi 3, although that might have been been my fritzy logic analyzer; I've since had no problems generating servo style pulses on all four pins with libbcm2835 on the Pi 3, which are only 50 Hz but rely on microsecond scale differences in width (I use a divisor of 8 and a range of 48000 for that).